Hm. Just a thought. Suppose we have a super-fast algorithm of finding intersections in 2D. Then, for 3D case, we can 1) project everything onto a plane — any one, XOY for example; 2) check for intersections on that plane; 3) then in 3D, check for intersections only for those pairs of edges which are intersecting in the projection.

How new implementation handle such coner cases as overlapping two edges when point of one edges lays on the edges of another and overlapping of two collinear edges?

Hm. Just a thought. Suppose we have a super-fast algorithm of finding intersections in 2D. Then, for 3D case, we can 1) project everything onto a plane — any one, XOY for example; 2) check for intersections on that plane; 3) then in 3D, check for intersections only for those pairs of edges which are intersecting in the projection.

That could, there would be more corner cases to check..... if two segments lay on a perpendicular plane to the projected plane the may look that they overlap but they could be intersecting...

The complicated part of this algorithms are the corner cases so I'm not sure it would be so simple to implement.

This could be interesting if someone finds the peace to implement it.
https://www.geometrictools.com/Documentation/DistanceLine3Line3.pdf

I thinks the best solution would be a numba implementation (numba is pip-installable) or a gpu one... in a perfect world

How new implementation handle such coner cases as overlapping two edges when point of one edges lays on the edges of another and overlapping of two collinear edges?

It does not handle them... :( at this point the remove doubles toggle is needed to remove the double point that it is created on the only-touching corner case... the overlapping edges are ignored.
I think I'm going to leave it as another option, not removing the old mathutils implementation, everything will be written on the documentation

Another option could be to ignore the only-touching case, that would not produce double points...

Maybe a extra toggle to let the user chose could be added...

Ok, finally the new algorithm was added as a new option of the 2d algorithms.
The previous Alg_1 have been optimised to be faster (it had a imperfect double removal that was slowing the algorithm to give a result that had some doubled vertices) The variation of this algorithm took this node to the mk3 version
Now Alg_1 and Np algorithms have a toggle to include On touch corner case (at the cost of generating doubled vertices that have to be removed with the r-doubles switch)

Also a new 3d intersection algorithm (in numpy) replaces the old one.
In this case the old one was taking 17 secs the new one 0.1 secs (170x)